#include "fit_pro.hh"

using namespace std;

fit_pro::fit_pro(string input_filename, string align_filename, string rotate_filename, string range_filename, float fit_rmse, bool is_have_target_det, int target_det_idx, int use_set_layer_used, double eff_sigma, double ADC_LSB, string save_path, bool is_write_waitalign_data, bool is_write_event_json, bool is_fit_control_angle, double angle_deg, string save_path_event, bool is_region_eff, int region_number)
{
    this->fit_rmse = fit_rmse;
    this->is_have_target_det = is_have_target_det;
    this->target_det_idx = target_det_idx;
    this->use_set_layer_used = use_set_layer_used;
    this->eff_sigma = eff_sigma;
    this->ADC_LSB = ADC_LSB;
    this->save_path = save_path;
    this->save_path_event = save_path_event;
    this->is_write_waitalign_data = is_write_waitalign_data;
    this->is_write_event_json = is_write_event_json;
    this->is_fit_control_angle = is_fit_control_angle;
    this->angle_deg = angle_deg;
    this->is_region_eff = is_region_eff;
    this->region_number = region_number;

    target_region_effx.resize(sqrt(region_number));
    target_region_effy.resize(sqrt(region_number));
    target_region_par.resize(region_number);
    ref_region_par.resize(region_number);
    region_par_amp.resize(region_number);
    ref_region_effx.resize(sqrt(region_number));
    ref_region_effy.resize(sqrt(region_number));

    // 前提条件
    if (!read_dec_file(align_filename, alignment))
    {
        return;
    }
    if (!read_dec_file(range_filename, fit_range))
    {
        return;
    }
    if (!read_dec_file(rotate_filename, rotate_data))
    {
        return;
    }

    // 读取文件，拟合运行
    int data_number = read_root_file(input_filename);
    is_file_get = data_number > 0;
    if (is_file_get)
    {
        set_layer_used();
        read_alignment();
        read_rotate();
        run();
    }
    else
    {
        cout << cRED << "Data in file initial failed!" << cRESET << endl;
        return;
    }
}

fit_pro::~fit_pro()
{
}

int fit_pro::read_root_file(string data_infilename)
{
    hit_data_file = new TFile(data_infilename.c_str());
    if (hit_data_file->IsZombie() || !hit_data_file->IsOpen())
    {
        cout << cRED << "[ERROR]: File" << data_infilename << " open failed. Please check the file" << cRESET << endl;
        return -1;
    }
    hit_data_tree = (TTree *)hit_data_file->Get("fec_hit_data");
    if (hit_data_tree == nullptr)
    {
        cout << cRED << "[ERROR]: Cannot get TTree* fec_hit_data in file " << data_infilename << cRESET << endl;
        return -1;
    }
    hit_data_tree->SetBranchAddress("trigger_id_all", &trigger_id, &b_triggera_id);
    hit_data_tree->SetBranchAddress("total_layer_num", &total_layer_num, &b_total_layer_num);
    hit_data_tree->GetEntry(0);
    valid_det_num = total_layer_num;
    for (int i = 0; i < valid_det_num; i++)
    {

        string dec_name = "dector" + to_string(i);
        hit_data_tree->SetBranchAddress(dec_name.c_str(), &detector_data[i].x_nhits, &b_detector_data[i]);
    }
    cout << cGREEN << "TTree in TFile " << data_infilename << " get!" << cRESET << endl;
    total_entried_num = hit_data_tree->GetEntries();
    if (total_entried_num != 0)
    {
        hit_data_tree->GetEntry(0);
    }
    return total_entried_num;
}

bool fit_pro::run()
{
    fit_vld_numx = 0;
    fit_vld_numy = 0;
    target_vld_numx = 0;
    target_vld_numy = 0;
    // hit_data_tree->GetEntry(0);

    for (int hidx = 0; hidx < total_entried_num; hidx++)
    {
        hit_data_tree->GetEntry(hidx);
        bool is_un_zero = false;
        for (int i = 0; i < valid_det_num; i++)
        {
            if (detector_data[i].z > 0)
            {
                is_un_zero = true;
            }
        }
        if (is_un_zero)
        {
            for (int i = 0; i < valid_det_num; i++)
            {
                if (!user_layer[i])
                {
                    continue;
                }
                if (is_have_target_det)
                {
                    if (i != target_det_idx)
                    {
                        dec_idx.push_back(i);
                        dec_construct.push_back(detector_data[i].z);
                    }
                    else
                    {
                        fit_location = detector_data[i].z;
                    }
                }
                else
                {
                    z[i] = detector_data[i].z;
                    dec_idx.push_back(i);
                    dec_construct.push_back(detector_data[i].z);
                }
            }
            break;
        }
    }

    // for (int i = 0; i < dec_construct.size(); i++)
    // {
    //     cout << dec_construct[i] << " ";
    // }
    // cout << endl;
    // cout << fit_location << endl;

    // 参考径迹探测器的数量
    //////////////////////////////
    det_layer_used = dec_idx.size();
    /////////////////////////////

    // 输出文件初始化
    /////////////////////////////
    if (!is_have_target_det)
    {
        write_root_file();
        valid_num = 0;
        event_id = 0;
        x_slope = 0;
        x_Intercept = 0;
        y_slope = 0;
        y_Intercept = 0;

        for (int i = 0; i < det_layer_used; i++)
        {
            detector_vld_data[i] = track_data();
        }
    }
    else if (is_have_target_det)
    {
        write_target_root_file();
        event_id = 0;
        valid_num = 0;
        det_idx = 0;

        x_slope = 0;
        x_Intercept = 0;
        x_rmse = 0;

        y_slope = 0;
        y_Intercept = 0;
        y_rmse = 0;

        x_pos = 0;
        y_pos = 0;
        z_pos = 0;

        x_amp = 0;
        y_amp = 0;

        x_pos_err = 0;
        y_pos_err = 0;

        x_strip_number = 0;
        y_strip_number = 0;
        is_x_meet = 0;
        is_y_meet = 0;
    }

    string filename = save_path + "fit_wait_aln.txt";
    ofstream out_wait_aln(filename, ios::app);
    string filename1 = save_path + "fit_valid_num.txt";
    ofstream num_vld(filename1, ios::app);
    /////////////////////////

    int one_percent = total_entried_num / rand() % 20 + 80;
    one_percent = (one_percent == 0) ? 1 : one_percent; // In case of total_entried_num < 100
    setbuf(stdout, NULL);

    cout << endl;

    cout << cRED << "[INFO]: total_entried_num: " << total_entried_num << cRESET << endl;
    bool is_write_title = true;
    for (int hit_idx = 0; hit_idx < total_entried_num; hit_idx++)
    {
        if (hit_idx % one_percent == 0)
        {
            cout << cCLR_LINE;
            double percent = (double)hit_idx / total_entried_num * 100.0;
            // Plot a progress bar
            cout << "[";
            int pos = 50 * percent / 100;
            for (int i = 0; i < 50; ++i)
            {
                if (i < pos)
                    cout << "=";
                else if (i == pos)
                    cout << ">";
                else
                    cout << " ";
            }
            // cout with 2 decimal points
            cout << fixed << setprecision(2);
            cout << "Processing " << percent << "%";
            // Restore the default setting
            cout << setprecision(6);
        }
        if (hit_idx == total_entried_num - 1)
        {
            cout << cCLR_LINE;
            cout << endl;
        }
        hit_data_tree->GetEntry(hit_idx);
        float *ref_det_chnx = new float[1000], *ref_det_chny = new float[1000], *target_det_chnx = new float[1000], *target_det_chny = new float[1000];
        float *target_det_ampx = new float[1000], *target_det_ampy = new float[1000];
        int *ref_det_dimx = new int[1000], *ref_det_dimy = new int[1000], *target_det_dimx = new int[1000], *target_det_dimy = new int[1000];
        float *ref_det_clux = new float[1000], *ref_det_cluy = new float[1000], *target_det_clux = new float[1000], *target_det_cluy = new float[1000];

        int xnum = 0;
        int ynum = 0;
        int ref_det_x_idx = 0, ref_det_y_idx = 0;
        int target_det_x_idx = 0, target_det_y_idx = 0;

        for (int i = 0; i < valid_det_num; i++)
        {
            if (!user_layer[i])
            {
                continue;
            }

            if (detector_data[i].x_nhits > 0)
            {
                if (is_have_target_det)
                {
                    if (i != target_det_idx)
                    {
                        xnum++;
                    }
                }
                else
                {
                    xnum++;
                }
            }

            if (detector_data[i].y_nhits > 0)
            {
                if (is_have_target_det)
                {
                    if (i != target_det_idx)
                    {
                        ynum++;
                    }
                }
                else
                {
                    ynum++;
                }
            }
        }

        if (xnum != det_layer_used || ynum != det_layer_used)
        {
            delete[] ref_det_chnx;
            delete[] ref_det_chny;
            delete[] target_det_chnx;
            delete[] target_det_chny;
            delete[] target_det_ampx;
            delete[] target_det_ampy;
            delete[] ref_det_dimx;
            delete[] ref_det_dimy;
            delete[] target_det_dimx;
            delete[] target_det_dimy;
            delete[] ref_det_clux;
            delete[] ref_det_cluy;
            delete[] target_det_clux;
            delete[] target_det_cluy;
            continue;
        }
        else
        {
            for (int i = 0; i < valid_det_num; i++)
            {
                if (!user_layer[i])
                {
                    continue;
                }

                if (detector_data[i].x_nhits > 0)
                {
                    for (int j = 0; j < detector_data[i].x_nhits; j++)
                    {
                        if (i == target_det_idx && is_have_target_det)
                        {
                            if (i == 9)
                            {
                                target_det_chnx[target_det_x_idx] = (detector_data[i].x[j] - 256) * cSTRIP_WIDTH;
                                target_det_ampx[target_det_x_idx] = abs(detector_data[i].x_amp[j]) * ADC_LSB;
                                target_det_clux[target_det_x_idx] = detector_data[i].x_chn_num[j];
                                target_det_dimx[target_det_x_idx] = i;
                                target_det_x_idx++;
                            }
                            else if (i == 5)
                            {
                                target_det_chnx[target_det_x_idx] = (detector_data[i].x[j] + 256) * cSTRIP_WIDTH;
                                target_det_ampx[target_det_x_idx] = abs(detector_data[i].x_amp[j]) * ADC_LSB;
                                target_det_clux[target_det_x_idx] = detector_data[i].x_chn_num[j];
                                target_det_dimx[target_det_x_idx] = i;
                                target_det_x_idx++;
                            }
                            else
                            {
                                target_det_chnx[target_det_x_idx] = (detector_data[i].x[j]) * cSTRIP_WIDTH;
                                target_det_ampx[target_det_x_idx] = abs(detector_data[i].x_amp[j]) * ADC_LSB;
                                target_det_clux[target_det_x_idx] = detector_data[i].x_chn_num[j];
                                target_det_dimx[target_det_x_idx] = i;
                                target_det_x_idx++;
                            }
                        }
                        else
                        {
                            if (i == 9)
                            {
                                ref_det_chnx[ref_det_x_idx] = (detector_data[i].x[j] - 256) * cSTRIP_WIDTH;
                                ref_det_clux[ref_det_x_idx] = detector_data[i].x_chn_num[j];
                                ref_det_dimx[ref_det_x_idx] = i;
                                ref_det_x_idx++;
                            }
                            else if (i == 5)
                            {
                                ref_det_chnx[ref_det_x_idx] = (detector_data[i].x[j] + 256) * cSTRIP_WIDTH;
                                ref_det_clux[ref_det_x_idx] = detector_data[i].x_chn_num[j];
                                ref_det_dimx[ref_det_x_idx] = i;
                                ref_det_x_idx++;
                            }
                            else
                            {
                                ref_det_chnx[ref_det_x_idx] = detector_data[i].x[j] * cSTRIP_WIDTH;
                                ref_det_clux[ref_det_x_idx] = detector_data[i].x_chn_num[j];
                                ref_det_dimx[ref_det_x_idx] = i;
                                ref_det_x_idx++;
                            }
                        }
                    }
                }

                if (detector_data[i].y_nhits > 0)
                {
                    for (int j = 0; j < detector_data[i].y_nhits; j++)
                    {
                        if (i == target_det_idx && is_have_target_det)
                        {
                            if (i == 9)
                            {
                                target_det_chny[target_det_y_idx] = (detector_data[i].y[j] - 256) * cSTRIP_WIDTH - 110;
                                target_det_ampy[target_det_y_idx] = abs(detector_data[i].y_amp[j]) * ADC_LSB;
                                target_det_cluy[target_det_y_idx] = detector_data[i].y_chn_num[j];
                                target_det_dimy[target_det_y_idx] = i;
                                target_det_y_idx++;
                            }
                            else if (i == 5)
                            {
                                target_det_chny[target_det_y_idx] = (detector_data[i].y[j] + 256) * cSTRIP_WIDTH;
                                target_det_ampy[target_det_y_idx] = abs(detector_data[i].y_amp[j]) * ADC_LSB;
                                target_det_cluy[target_det_y_idx] = detector_data[i].y_chn_num[j];
                                target_det_dimy[target_det_y_idx] = i;
                                target_det_y_idx++;
                            }
                            else
                            {
                                target_det_chny[target_det_y_idx] = (detector_data[i].y[j]) * cSTRIP_WIDTH;
                                target_det_ampy[target_det_y_idx] = abs(detector_data[i].y_amp[j]) * ADC_LSB;
                                target_det_cluy[target_det_y_idx] = detector_data[i].y_chn_num[j];
                                target_det_dimy[target_det_y_idx] = i;
                                target_det_y_idx++;
                            }
                        }
                        else
                        {
                            if (i == 9)
                            {
                                ref_det_chny[ref_det_y_idx] = (detector_data[i].y[j] - 256) * cSTRIP_WIDTH - 110;
                                ref_det_cluy[ref_det_y_idx] = detector_data[i].y_chn_num[j];
                                ref_det_dimy[ref_det_y_idx] = i;
                                ref_det_y_idx++;
                            }
                            else if (i == 5)
                            {
                                ref_det_chny[ref_det_y_idx] = (detector_data[i].y[j] + 256) * cSTRIP_WIDTH;
                                ref_det_cluy[ref_det_y_idx] = detector_data[i].y_chn_num[j];
                                ref_det_dimy[ref_det_y_idx] = i;
                                ref_det_y_idx++;
                            }
                            else
                            {
                                ref_det_chny[ref_det_y_idx] = detector_data[i].y[j] * cSTRIP_WIDTH;
                                ref_det_cluy[ref_det_y_idx] = detector_data[i].y_chn_num[j];
                                ref_det_dimy[ref_det_y_idx] = i;
                                ref_det_y_idx++;
                            }
                        }
                    }
                }
            }
        }

        // 数据组合筛选过滤
        int wait_fit_num_x = 0, wait_fit_num_y = 0;
        float *wait_fit_data_rx = nullptr, *wait_fit_data_ry = nullptr, *wait_fit_data_cluster_rx = nullptr, *wait_fit_data_cluster_ry = nullptr;

        bool x = fit_data_select(ref_det_chnx, ref_det_clux, ref_det_dimx, ref_det_x_idx, wait_fit_num_x, wait_fit_data_rx, wait_fit_data_cluster_rx);
        bool y = fit_data_select(ref_det_chny, ref_det_cluy, ref_det_dimy, ref_det_y_idx, wait_fit_num_y, wait_fit_data_ry, wait_fit_data_cluster_ry);

        if (x && y)
        {
        }
        else
        {
            delete[] ref_det_chnx;
            delete[] ref_det_chny;
            delete[] ref_det_dimx;
            delete[] ref_det_dimy;
            delete[] ref_det_clux;
            delete[] ref_det_cluy;
            delete[] target_det_chnx;
            delete[] target_det_chny;
            delete[] target_det_ampx;
            delete[] target_det_ampy;
            delete[] target_det_dimx;
            delete[] target_det_dimy;
            delete[] target_det_clux;
            delete[] target_det_cluy;
            continue;
        }

        // 旋转数据 等待拟合
        float *wait_fit_data_ox = nullptr, *wait_fit_data_oy = nullptr, *wait_fit_data_oz = nullptr, *wait_fit_data_cluster_ox = nullptr, *wait_fit_data_cluster_oy = nullptr;
        int wait_fit_num_ox = 0, wait_fit_num_oy = 0, all_size = 0;
        if (is_rotate)
        {
            data_rotate(wait_fit_num_x, wait_fit_num_y, wait_fit_data_rx, wait_fit_data_ry, wait_fit_data_cluster_rx, wait_fit_data_cluster_ry, wait_fit_data_ox, wait_fit_data_oy, wait_fit_data_oz, wait_fit_data_cluster_ox, wait_fit_data_cluster_oy, all_size);
            wait_fit_num_ox = all_size / det_layer_used;
            wait_fit_num_oy = all_size / det_layer_used;
        }
        else
        {
            wait_fit_num_ox = wait_fit_num_x;
            wait_fit_num_oy = wait_fit_num_y;
            wait_fit_data_ox = wait_fit_data_rx;
            wait_fit_data_oy = wait_fit_data_ry;
            wait_fit_data_oz = new float[det_layer_used];
            for (int i = 0; i < det_layer_used; i++)
            {
                wait_fit_data_oz[i] = dec_construct[i];
            }
            wait_fit_data_cluster_ox = wait_fit_data_cluster_rx;
            wait_fit_data_cluster_oy = wait_fit_data_cluster_ry;
        }

        // 拟合数据
        float fit_hit_x = 0, fit_hit_y = 0;
        float fit_slopex = 0, fit_slopey = 0;
        float fit_interceptx = 0, fit_intercepty = 0;

        if (is_rotate)
        {
            det_data_fitv2(wait_fit_num_ox, wait_fit_data_ox, wait_fit_data_oy, wait_fit_data_oz, wait_fit_data_cluster_ox, wait_fit_data_cluster_oy, fit_hit_x, fit_hit_y, fit_slopex, fit_slopey, fit_interceptx, fit_intercepty);
        }
        else
        {
            det_data_fit(wait_fit_num_ox, wait_fit_data_ox, wait_fit_data_oz, wait_fit_data_cluster_ox, true, fit_hit_x, fit_slopex, fit_interceptx);
            det_data_fit(wait_fit_num_oy, wait_fit_data_oy, wait_fit_data_oz, wait_fit_data_cluster_oy, false, fit_hit_y, fit_slopey, fit_intercepty);
        }

        if (!is_have_target_det && fec_fit_datax.size() != 0 && fec_fit_datay.size() != 0)
        {
            vector<float> temp_data;
            if (fec_fit_datax[0][0] == fec_fit_datay[0][0])
            {
                event_id = fec_fit_datax[0][0];
                x_slope = fec_fit_datax[0][1];
                x_Intercept = fec_fit_datax[0][2];
                y_slope = fec_fit_datay[0][1];
                y_Intercept = fec_fit_datay[0][2];

                for (int i = 0; i < det_layer_used; i++)
                {
                    detector_vld_data[i].x = fec_fit_datax[0][i + 3];
                    detector_vld_data[i].y = fec_fit_datay[0][i + 3];
                    detector_vld_data[i].z = fec_fit_datax[0][i + 3 + det_layer_used];
                }

                adas_track_tree->Fill();

                if (is_write_event_json)
                {
                    trig_event = fec_fit_datax[0][0];
                    kx = fec_fit_datax[0][1];
                    ky = fec_fit_datay[0][1];
                    bx = fec_fit_datax[0][2];
                    by = fec_fit_datay[0][2];
                    for (int i = 0; i < det_layer_used; i++)
                    {
                        x_strip[i] = fec_fit_datax[0][i + 3];
                        y_strip[i] = fec_fit_datay[0][i + 3];
                        x_strip_num[i] = fec_fit_datax[0][i + 3 + 2 * det_layer_used];
                        y_strip_num[i] = fec_fit_datay[0][i + 3 + 2 * det_layer_used];
                    }
                    write_event_json_file(event_id);
                    write_status_json_file();
                }

                if (is_write_waitalign_data)
                {
                    out_wait_aln << event_id << " ";
                    for (int i = 0; i < det_layer_used; i++)
                    {
                        out_wait_aln << detector_vld_data[i].x << " ";
                        out_wait_aln << detector_vld_data[i].y << " ";
                        out_wait_aln << detector_vld_data[i].z << " ";
                    }
                    out_wait_aln << endl;
                }
                valid_num++;
            }
            fec_fit_datax.clear();
            fec_fit_datay.clear();
        }

        if (1)
        {
            delete[] wait_fit_data_ox;
            delete[] wait_fit_data_oy;
            delete[] wait_fit_data_oz;
            delete[] wait_fit_data_cluster_ox;
            delete[] wait_fit_data_cluster_oy;
        }
        float target_ampx = 0, target_chnx = 0, target_clux = 0, target_x_err = 0;
        float target_ampy = 0, target_chny = 0, target_cluy = 0, target_y_err = 0;
        int target_region_idx_x = 0, target_region_idx_y = 0, par_region_idx = 0;

        if (is_region_eff)
        {
            int n = sqrt(region_number);
            if (n * n != region_number)
            {
                cout << cRED << "[ERROR]: region_number must be a perfect square!" << cRESET << endl;
                return false;
            }

            double region_size = (fit_range[1] - fit_range[0]) / n;

            if (fit_hit_x != 0 && fit_hit_y != 0)
            {
                int region_idx_x = (int)((fit_hit_x - fit_range[0]) / region_size);
                target_region_idx_x = region_idx_x;
                ref_region_effx[region_idx_x].push_back(trigger_id);

                int region_idx_y = (int)((fit_hit_y - fit_range[0]) / region_size);
                target_region_idx_y = region_idx_y;
                ref_region_effy[region_idx_y].push_back(trigger_id);
                par_region_idx = (region_idx_x)*n + region_idx_y;
                ref_region_par[par_region_idx].push_back(trigger_id);
            }
        }

        if (fit_hit_x != 0 && is_have_target_det && fit_hit_y != 0)
        {
            float *t_r_x = nullptr, *t_r_y = nullptr, *t_r_z = nullptr, *t_c_x = nullptr, *t_c_y = nullptr, *t_a_x = nullptr, *t_a_y = nullptr, *t_r_zx = nullptr, *t_r_zy = nullptr;
            int t_numx = 0, t_numy = 0, t_num = 0;

            fit_vld_numx++;
            fit_vld_numy++;
            if (target_det_x_idx > 0 && target_det_y_idx > 0)
            {

                if (is_rotate)
                {
                    target_rotate(target_det_x_idx, target_det_y_idx, target_det_chnx, target_det_chny, target_det_clux, target_det_cluy, target_det_ampx, target_det_ampy, t_r_x, t_r_y, t_r_z, t_c_x, t_c_y, t_a_x, t_a_y, t_num);
                    t_numx = t_num;
                    t_numy = t_num;
                    t_r_zx = new float[1000];
                    t_r_zy = new float[1000];
                    for (int i = 0; i < t_num; i++)
                    {
                        float temp = t_r_z[i];
                        t_r_zx[i] = temp;
                        t_r_zy[i] = temp;
                    }
                }
                else
                {
                    t_r_zx = new float[1000];
                    t_r_zy = new float[1000];
                    for (int i = 0; i < target_det_x_idx; i++)
                    {
                        t_r_zx[i] = fit_location;
                    }
                    for (int i = 0; i < target_det_y_idx; i++)
                    {
                        t_r_zy[i] = fit_location;
                    }

                    t_numx = target_det_x_idx;
                    t_numy = target_det_y_idx;
                    t_r_x = target_det_chnx;
                    t_r_y = target_det_chny;
                    t_c_x = target_det_clux;
                    t_c_y = target_det_cluy;
                    t_a_x = target_det_ampx;
                    t_a_y = target_det_ampy;
                }

                calculateSpatial(t_r_zx, t_r_x, t_a_x, t_c_x, t_numx, fit_hit_x, fit_slopex, fit_interceptx, tx, target_ampx, target_chnx, target_clux, target_x_err, target_y_err, true);
                calculateSpatial(t_r_zy, t_r_y, t_a_y, t_c_y, t_numy, fit_hit_y, fit_slopey, fit_intercepty, ty, target_ampy, target_chny, target_cluy, target_x_err, target_y_err, false);

                if (target_ampx != 0)
                {
                    target_vld_numx++;
                    if (is_region_eff)
                    {
                        target_region_effx[target_region_idx_x].push_back(target_ampx);
                    }
                }

                if (target_ampy != 0)
                {
                    target_vld_numy++;
                    if (is_region_eff)
                    {
                        target_region_effy[target_region_idx_y].push_back(target_ampy);
                    }
                }
            }
            else if (!((target_det_x_idx > 0) && (target_det_y_idx > 0)))
            {

                if (is_rotate)
                {
                    if (target_det_x_idx > 0)
                    {
                        float *temp_chny = new float[1000];
                        float *temp_ampy = new float[1000];
                        float *temp_cluy = new float[1000];
                        temp_chny[0] = fit_hit_y;
                        temp_ampy[0] = 0;
                        temp_cluy[0] = 0;
                        target_rotate(target_det_x_idx, 1, target_det_chnx, temp_chny, target_det_clux, temp_cluy, target_det_ampx, temp_ampy, t_r_x, t_r_y, t_r_z, t_c_x, t_c_y, t_a_x, t_a_y, t_num);
                        t_numx = t_num;
                        t_numy = 0;
                        t_r_zx = new float[1000];
                        t_r_zy = new float[1000];
                        for (int i = 0; i < t_num; i++)
                        {
                            float temp = t_r_z[i];
                            t_r_zx[i] = temp;
                            t_r_zy[i] = temp;
                        }
                        // cout << endl;
                        // cout << "fit_hit_y: " << fit_hit_y << endl;
                        // for (int i = 0; i < t_num; i++)
                        // {
                        //     cout << t_r_zx[i] << " ";
                        //     cout << t_r_x[i] << " ";
                        //     cout << t_r_y[i] << " ";
                        // }
                        // cout << endl;
                    }
                    else if (target_det_y_idx > 0)
                    {
                        float *temp_chnx = new float[1000];
                        float *temp_ampx = new float[1000];
                        float *temp_clux = new float[1000];
                        temp_chnx[0] = fit_hit_x;
                        temp_ampx[0] = 0;
                        temp_clux[0] = 0;
                        target_rotate(1, target_det_y_idx, temp_chnx, target_det_chny, temp_clux, target_det_cluy, temp_ampx, target_det_ampy, t_r_x, t_r_y, t_r_z, t_c_x, t_c_y, t_a_x, t_a_y, t_num);
                        t_numx = 0;
                        t_numy = t_num;
                        t_r_zx = new float[1000];
                        t_r_zy = new float[1000];
                        for (int i = 0; i < t_num; i++)
                        {
                            float temp = t_r_z[i];
                            t_r_zx[i] = temp;
                            t_r_zy[i] = temp;
                        }
                    }
                }
                else
                {
                    t_r_zx = new float[1000];
                    t_r_zy = new float[1000];
                    for (int i = 0; i < target_det_x_idx; i++)
                    {
                        t_r_zx[i] = fit_location;
                    }
                    for (int i = 0; i < target_det_y_idx; i++)
                    {
                        t_r_zy[i] = fit_location;
                    }

                    t_numx = target_det_x_idx;
                    t_numy = target_det_y_idx;
                    t_r_x = target_det_chnx;
                    t_r_y = target_det_chny;
                    t_c_x = target_det_clux;
                    t_c_y = target_det_cluy;
                    t_a_x = target_det_ampx;
                    t_a_y = target_det_ampy;
                }

                if (t_numx > 0)
                {
                    calculateSpatial(t_r_zx, t_r_x, t_a_x, t_c_x, t_numx, fit_hit_x, fit_slopex, fit_interceptx, tx, target_ampx, target_chnx, target_clux, target_x_err, target_y_err, true);
                }
                if (t_numy > 0)
                {
                    calculateSpatial(t_r_zy, t_r_y, t_a_y, t_c_y, t_numy, fit_hit_y, fit_slopey, fit_intercepty, ty, target_ampy, target_chny, target_cluy, target_x_err, target_y_err, false);
                }

                if (target_ampx != 0)
                {
                    target_vld_numx++;
                    if (is_region_eff)
                    {
                        target_region_effx[target_region_idx_x].push_back(target_ampx);
                    }
                }

                if (target_ampy != 0)
                {
                    target_vld_numy++;
                    if (is_region_eff)
                    {
                        target_region_effy[target_region_idx_y].push_back(target_ampy);
                    }
                }
            }

            delete[] t_r_x;
            delete[] t_r_y;
            delete[] t_r_z;
            delete[] t_c_x;
            delete[] t_c_y;
            delete[] t_a_x;
            delete[] t_a_y;
            delete[] t_r_zx;
            delete[] t_r_zy;
        }

        if (target_ampx != 0 && target_ampy != 0 && is_have_target_det)
        {
            ampall.push_back(target_ampx + target_ampy);
            x_all.push_back(target_chnx);
            y_all.push_back(target_chny);
            target_region_par[par_region_idx].push_back(target_ampx);
            region_par_amp[par_region_idx].push_back(target_ampx + target_ampy);
        }

        if (fit_hit_x != 0 && is_have_target_det && fit_hit_y != 0)
        {

            event_id = trigger_id;
            layer_num = total_layer_num;
            det_idx = target_det_idx + 1;
            x_slope = fec_fit_datax[0][1];
            x_Intercept = fec_fit_datax[0][2];
            x_rmse = fec_fit_datax[0][3];
            y_slope = fec_fit_datay[0][1];
            y_Intercept = fec_fit_datay[0][2];
            y_rmse = fec_fit_datay[0][3];

            is_x_meet = (target_ampx != 0);
            is_y_meet = (target_ampy != 0);
            z_pos = fit_location + tz;
            if (is_x_meet)
            {
                x_pos = target_chnx;
                x_amp = target_ampx;
                x_strip_number = target_clux;
                x_pos_err = target_x_err;
            }
            else
            {
                x_pos = 0;
                x_amp = 0;
                x_strip_number = 0;
                x_pos_err = 0;
            }

            if (is_y_meet)
            {
                y_pos = target_chny;
                y_amp = target_ampy;
                y_strip_number = target_cluy;
                y_pos_err = target_y_err;
            }
            else
            {
                y_pos = 0;
                y_amp = 0;
                y_strip_number = 0;
                y_pos_err = 0;
            }

            adas_target_tree->Fill();
            if (is_write_waitalign_data)
            {
                ofstream out_fit_data;
                string fit_data_filename = save_path + "target_layer" + to_string(target_det_idx + 1) + "_align_fit_data.txt";
                out_fit_data.open(fit_data_filename, ios::app);
                if (out_fit_data.is_open())
                {
                    if (is_write_title)
                    {
                        out_fit_data << "trigger_id: " << "x_slope: " << "x_Intercept: " << "x_rmse: " << "y_slope: " << "y_Intercept: " << "y_rmse: " << "hit_idx: " << endl;
                        is_write_title = false;
                    }
                    for (size_t i = 0; i < fec_fit_datax.size(); i++)
                    {
                        out_fit_data << trigger_id << " ";
                        for (size_t j = 1; j < fec_fit_datax[i].size(); j++)
                        {
                            out_fit_data << fec_fit_datax[i][j] << " ";
                        }
                        for (size_t j = 1; j < fec_fit_datay[i].size(); j++)
                        {
                            out_fit_data << fec_fit_datay[i][j] << " ";
                        }
                        out_fit_data << hit_idx << " ";
                        out_fit_data << endl;
                    }
                    out_fit_data.close();
                }
            }
        }
        fec_fit_datax.clear();
        fec_fit_datay.clear();
    }

    if (is_have_target_det)
    {
        if (is_write_waitalign_data)
        {
            string feename;
            if (ADC_LSB == 1)
            {
                feename = "STAGE_";
            }
            else
            {
                feename = "ADAS1128_";
            }
            string x_name = feename + to_string(target_det_idx + 1) + "_x.txt";

            ofstream outallx(save_path + x_name, ios::app);
            for (int i = 0; i < ampx.size(); i++)
            {
                // outallx <<"amp: "<< ampx[i]/28<< " rms: "<<rmsx[i]<<" chnx: "<<chnx[i] << endl;
                outallx << targe_vld_idx[i] << " " << ampx[i] << " " << rmsx[i] << " " << chnx[i] << endl;
            }
            outallx.close();

            string y_name = feename + to_string(target_det_idx + 1) + "_y.txt";

            ofstream outally(save_path + y_name, ios::app);
            for (int i = 0; i < ampy.size(); i++)
            {
                // outally <<"amp: "<< ampy[i]/28<<" rms: "<<rmsy[i]<<" chny: "<<chny[i] << endl;
                outally << targe_vld_idy[i] << " " << ampy[i] << " " << rmsy[i] << " " << chny[i] << endl;
            }
            outally.close();
        }
        out_result();
        draw_spatial(rmsx, true);
        draw_spatial(rmsy, false);
        draw_energy(ampx, true);
        draw_energy(ampy, false);
        if (!is_fit_control_angle)
        {
            draw_cluster(clusterx, true);
            draw_cluster(clustery, false);
            draw_all_energy(ampall);
            draw_hit_range_xy();
            draw_hit_strip(x_all, true);
            draw_hit_strip(y_all, false);
        }
        if (is_region_eff)
        {
            draw_region_efficiency();
            // draw_region_amp();
        }
        cout << target_vld_numx << " " << target_vld_numy << endl;
        cout << fit_vld_numx << " " << fit_vld_numy << endl;

        string effname = "x_eff_" + to_string(target_det_idx + 1);
        float effx = (float)target_vld_numx / (float)fit_vld_numx;
        outall << effname << " " << effx << endl;

        cout << cRED << "detector " << target_det_idx + 1 << " x direction efficiency: " << effx << cRESET << endl;

        effname = "y_eff_" + to_string(target_det_idx + 1);
        float effy = (float)target_vld_numy / (float)fit_vld_numy;
        outall << effname << " " << effy << endl;

        cout << cRED << "detector " << target_det_idx + 1 << " y direction efficiency: " << effy << cRESET << endl;
        complete_adding();
    }
    else
    {
        cout << cRED << "Total valid event number: " << valid_num << cRESET << endl;
        num_vld << valid_num << endl;
        if (!is_write_event_json)
        {
            cout << cBLUE << "det wait alignment data output sucess" << cRESET << endl;
        }
        else
        {
            cout << cBLUE << "event json write success" << save_path_event << cRESET << endl;
        }
        out_wait_aln.close();
        complete_adding();
    }

    cout << cBLUE << "Data all calculate success!" << cRESET << endl;
    cout << endl;
    return true;
}

bool fit_pro::fit_data_select(float *event_det_chn, float *event_det_cluster, int *event_det_idx, int event_det_size, int &assembly_num, float *&fit_data, float *&fit_cluster)
{
    // 1. 统计每个数字的出现次数
    std::map<int, int> counts;
    for (int i = 0; i < event_det_size; ++i)
    {
        counts[event_det_idx[i]]++;
    }

    // 3. 计算总组合数
    assembly_num = 1;
    for (int i = 0; i < det_layer_used; ++i)
    {
        assembly_num *= counts[dec_idx[i]];
    }

    // 4. 计算新数组的总长度
    int total_data_points = assembly_num * det_layer_used;

    // 5. 初始化偏移量数组
    std::vector<int> offsets(det_layer_used, 0);
    // 创建 data_map,  预先分配空间
    std::map<int, std::vector<int>> data_map;
    for (auto const &[num, count] : counts)
    {
        data_map[num].reserve(count);
    }

    for (int i = 0; i < event_det_size; ++i)
    {
        data_map[event_det_idx[i]].push_back(i);
    }
    int final_idx = 0;

    fit_data = new float[total_data_points];
    fit_cluster = new float[total_data_points];

    // 7. 生成所有组合
    for (int comb_idx = 0; comb_idx < assembly_num; ++comb_idx)
    {
        float *current_data = new float[det_layer_used];
        float *current_cluster = new float[det_layer_used];
        int current_data_idx = 0;

        for (int i = 0; i < det_layer_used; ++i)
        {
            int num = dec_idx[i];
            int data_index = -1;
            if (offsets[i] < data_map[num].size())
            {
                data_index = data_map[num][offsets[i]];
                current_data[current_data_idx] = event_det_chn[data_index];
                current_cluster[current_data_idx] = event_det_cluster[data_index];
                current_data_idx++;
            }
        }

        // 更新 offsets
        offsets[0]++;
        for (int i = 0; i < det_layer_used - 1; i++)
        {
            if (offsets[i] >= counts[dec_idx[i]])
            {
                offsets[i] = 0;
                offsets[i + 1]++;
            }
            else
            {
                break;
            }
        }

        int size = 0;
        float *temp_data = nullptr;
        float *temp_cluster = nullptr;

        temp_data = new float[det_layer_used];
        temp_cluster = new float[det_layer_used];

        if (!data_filter(1, current_data, current_cluster, size, temp_data, temp_cluster))
        {
            delete[] temp_data;
            delete[] temp_cluster;
            delete[] current_data;
            delete[] current_cluster;
            continue;
        }
        else
        {
            for (int i = 0; i < det_layer_used; i++)
            {
                fit_data[final_idx] = temp_data[i];
                fit_cluster[final_idx] = temp_cluster[i];
                final_idx++;
            }
            delete[] temp_data;
            delete[] temp_cluster;
            delete[] current_data;
            delete[] current_cluster;
        }
    }

    if (final_idx == 0)
    {
        return false;
    }
    else
    {
        assembly_num = final_idx / det_layer_used;
        return true;
    }
}

bool fit_pro::data_filter(int num, float *wait_fliter_data, float *wait_fliter_cluster, int &fliter_size, float *&fliter_over_data, float *&fliter_over_cluster)

{
    float rmse = 0;
    int data_size = det_layer_used;
    fliter_size = 0;
    int k = 0;

    float select_vth = 0;
    if (is_have_target_det)
    {
        select_vth = 30;
    }
    else
    {
        select_vth = fit_rmse;
    }

    fliter_over_data = new float[num * det_layer_used];
    fliter_over_cluster = new float[num * det_layer_used];
    for (int j = 0; j < num; j++)
    {
        float *current_fec_data_combine = wait_fliter_data + k;
        float *current_fec_data_combine_cluster = wait_fliter_cluster + k;

        k += data_size;
        VectorXd Y(data_size);
        MatrixXd X(data_size, 2);
        for (int n = 0; n < data_size; n++)
        {
            X(n, 0) = 1.0;
            X(n, 1) = dec_construct[n];
            Y(n) = current_fec_data_combine[n];
        }
        VectorXd beta = X.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Y);
        rmse = 0;
        for (int z = 0; z < data_size; z++)
        {
            float fit_value = beta(0) + beta(1) * X(z, 1);
            rmse += (Y(z) - fit_value) * (Y(z) - fit_value);
        }
        rmse = sqrt(rmse / data_size);
        // cout << rmse << endl;

        if (rmse < select_vth)
        {
            for (int z = 0; z < data_size; z++)
            {
                fliter_over_data[fliter_size] = current_fec_data_combine[z];
                fliter_over_cluster[fliter_size] = current_fec_data_combine_cluster[z];
                fliter_size++;
            }
        }
    }

    if (fliter_size == 0)
    {
        return false;
    }
    else
    {
        return true;
    }
}

void fit_pro::data_rotate(int numx, int numy, float *wait_rotatex, float *wait_rotatey, float *wait_clusterx, float *wait_clustery, float *&det_rotatex, float *&det_rotatey, float *&det_rotatez, float *&det_clusterx, float *&det_clustery, int &output_size)
{

    output_size = numx * numy * det_layer_used;

    // 1. Allocate memory for all arrays
    det_rotatex = new float[output_size];
    det_rotatey = new float[output_size];
    det_rotatez = new float[output_size];
    det_clusterx = new float[output_size];
    det_clustery = new float[output_size];

    memset(det_rotatex, 0, output_size * sizeof(float));
    memset(det_rotatey, 0, output_size * sizeof(float));
    memset(det_rotatez, 0, output_size * sizeof(float));
    memset(det_clusterx, 0, output_size * sizeof(float));
    memset(det_clustery, 0, output_size * sizeof(float));

    int mm = 0;
    // 使用指针的方式
    for (int i = 0; i < numx; i++)
    {
        for (int j = 0; j < numy; j++)
        {
            for (int k = 0; k < det_layer_used; k++)
            {
                float x = wait_rotatex[i * det_layer_used + k];
                float y = wait_rotatey[j * det_layer_used + k];
                float z = dec_construct[k];

                Matrix<float, 3, 1> detxyz;
                detxyz << x, y, z;

                Matrix<float, 3, 1> detxyz_rotated = rotate_matrix[k] * detxyz;

                det_rotatex[mm] = detxyz_rotated(0);
                det_rotatey[mm] = detxyz_rotated(1);
                det_rotatez[mm] = detxyz_rotated(2);
                det_clusterx[mm] = wait_clusterx[i * det_layer_used + k];
                det_clustery[mm] = wait_clustery[j * det_layer_used + k];

                mm++;
            }
        }
    }
}

void fit_pro::target_rotate(int numx, int numy, float *wait_rotatex, float *wait_rotatey, float *wait_clusterx, float *wait_clustery, float *wait_ampx, float *wait_ampy, float *&det_rotatex, float *&det_rotatey, float *&det_rotatez, float *&det_clusterx, float *&det_clustery, float *&det_ampx, float *&det_ampy, int &output_size)
{
    output_size = numx * numy;

    // 1. Allocate memory for all arrays
    det_rotatex = new float[output_size];
    det_rotatey = new float[output_size];
    det_rotatez = new float[output_size];
    det_clusterx = new float[output_size];
    det_clustery = new float[output_size];
    det_ampx = new float[output_size];
    det_ampy = new float[output_size];

    memset(det_rotatex, 0, output_size * sizeof(float));
    memset(det_rotatey, 0, output_size * sizeof(float));
    memset(det_rotatez, 0, output_size * sizeof(float));
    memset(det_clusterx, 0, output_size * sizeof(float));
    memset(det_clustery, 0, output_size * sizeof(float));
    memset(det_ampx, 0, output_size * sizeof(float));
    memset(det_ampy, 0, output_size * sizeof(float));

    int mm = 0;
    // 使用指针的方式
    for (int i = 0; i < numx; i++)
    {
        for (int j = 0; j < numy; j++)
        {

            float x = wait_rotatex[i];
            float y = wait_rotatey[j];
            float z = fit_location;

            Matrix<float, 3, 1> detxyz;
            detxyz << x, y, z;

            Matrix<float, 3, 1> detxyz_rotated = rotate_matrix_target * detxyz;

            det_rotatex[mm] = detxyz_rotated(0);
            det_rotatey[mm] = detxyz_rotated(1);
            det_rotatez[mm] = detxyz_rotated(2);
            det_clusterx[mm] = wait_clusterx[i];
            det_clustery[mm] = wait_clustery[j];
            det_ampx[mm] = wait_ampx[i];
            det_ampy[mm] = wait_ampy[j];
            mm++;
        }
    }
}

void fit_pro::det_data_fitv2(int num, float *fec_data_combine_x, float *fec_data_combine_y, float *fec_data_z, float *fec_data_combine_cluster_x, float *fec_data_combine_cluster_y, float &vld_fit_hit_x, float &vld_fit_hit_y, float &vld_fit_slope_x, float &vld_fit_slope_y, float &vld_fit_intercept_x, float &vld_fit_intercept_y)
{
    vld_fit_hit_x = 0;
    vld_fit_hit_y = 0;
    vld_fit_slope_x = 0;
    vld_fit_slope_y = 0;
    vld_fit_intercept_x = 0;
    vld_fit_intercept_y = 0;

    float fit_hit_x = 0;
    float fit_hit_y = 0;
    float slope_x = 0;
    float slope_y = 0;
    float Intercept_x = 0;
    float Intercept_y = 0;
    float rmse_x = 0;
    float rmse_y = 0;
    float rmse_all = 0;
    float rmse_all_min = fit_rmse;
    float fit_min_x = fit_rmse;
    float fit_min_y = fit_rmse;

    double tan_x = 0, tan_y = 0;
    double fit_angle_x = 0, fit_angle_y = 0;

    int data_size = det_layer_used;

    int k = 0;
    vector<float> out_data_x;
    vector<float> out_data_y;
    for (int j = 0; j < num; j++)
    {
        float *current_fec_data_combine_x = fec_data_combine_x + k;
        float *current_fec_data_combine_y = fec_data_combine_y + k;
        float *current_fec_data_z;
        if (is_rotate)
        {
            current_fec_data_z = fec_data_z + k;
        }
        else
        {
            current_fec_data_z = fec_data_z;
        }
        float *current_fec_data_combine_cluster_x = fec_data_combine_cluster_x + k;
        float *current_fec_data_combine_cluster_y = fec_data_combine_cluster_y + k;

        k += data_size;

        VectorXd Yx(data_size);
        VectorXd Yy(data_size);
        MatrixXd X(data_size, 2);

        for (int n = 0; n < data_size; n++)
        {
            X(n, 0) = 1.0;
            X(n, 1) = current_fec_data_z[n] + z_alignment[n];
            Yx(n) = current_fec_data_combine_x[n] + x_alignment[n];
            Yy(n) = current_fec_data_combine_y[n] + y_alignment[n];
        }
        VectorXd beta_x = X.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Yx);
        VectorXd beta_y = X.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Yy);

        // rmse_x
        rmse_x = 0;
        for (int z = 0; z < data_size; z++)
        {
            double fit_value_x = X.row(z) * beta_x; // 适用于多元
            rmse_x += (Yx(z) - fit_value_x) * (Yx(z) - fit_value_x);
        }
        rmse_x = std::sqrt(rmse_x / data_size);
        // rmse_y
        rmse_y = 0;
        for (int z = 0; z < data_size; z++)
        {
            double fit_value_y = X.row(z) * beta_y; // 适用于多元
            rmse_y += (Yy(z) - fit_value_y) * (Yy(z) - fit_value_y);
        }
        rmse_y = std::sqrt(rmse_y / data_size);
        if (is_fit_control_angle)
        {
            tan_x = beta_x(1);
            fit_angle_x = std::atan(std::abs(tan_x)) * (180 / M_PI);
            tan_y = beta_y(1);
            fit_angle_y = std::atan(std::abs(tan_y)) * (180 / M_PI);
        }
        else
        {
            fit_angle_x = 0;
            fit_angle_y = 0;
        }

        rmse_all = (rmse_x + rmse_y) / 2;

        if (rmse_x < fit_rmse && rmse_x < fit_min_x && fit_angle_x < angle_deg && rmse_y < fit_rmse && rmse_y < fit_min_y && fit_angle_y < angle_deg && rmse_all < fit_rmse && rmse_all < rmse_all_min)
        {
            fit_min_x = rmse_x;
            fit_min_y = rmse_y;
            slope_x = beta_x(1);
            slope_y = beta_y(1);
            Intercept_x = beta_x(0);
            Intercept_y = beta_y(0);
            rmse_all_min = rmse_all;

            for (int z = 0; z < data_size; z++)
            {
                out_data_x.push_back(Yx(z));
                out_data_y.push_back(Yy(z));
            }

            for (int z = 0; z < data_size; z++)
            {
                out_data_x.push_back(X(z, 1));
                out_data_y.push_back(X(z, 1));
            }

            if (is_write_event_json)
            {
                for (int z = 0; z < data_size; z++)
                {
                    out_data_x.push_back(current_fec_data_combine_cluster_x[z]);
                    out_data_y.push_back(current_fec_data_combine_cluster_y[z]);
                }
            }
            if (is_have_target_det)
            {
                fit_hit_x = (Intercept_x + slope_x * (fit_location + z_alignment[0]));
                fit_hit_y = (Intercept_y + slope_y * (fit_location + z_alignment[0]));
            }
            else
            {
                fit_hit_x = 0;
                fit_hit_y = 0;
            }
        }
    }

    if (fit_hit_x != 0 && fit_min_x < fit_rmse && fit_hit_x >= fit_range[0] && fit_hit_x <= fit_range[1] && is_have_target_det)
    {
        vld_fit_hit_x = fit_hit_x;
        vld_fit_slope_x = slope_x;
        vld_fit_intercept_x = Intercept_x;
        vector<float> temp_data = {(float)trigger_id, slope_x, Intercept_x, fit_min_x, fit_hit_x};
        // temp_data.insert(temp_data.end(), out_data.begin(), out_data.end());
        fec_fit_datax.push_back(temp_data);
    }
    else if (fit_min_x < fit_rmse && !is_have_target_det)
    {
        vector<float> temp_data = {(float)trigger_id, slope_x, Intercept_x};
        temp_data.insert(temp_data.end(), out_data_x.begin(), out_data_x.end());

        fec_fit_datax.push_back(temp_data);
    }

    if (fit_hit_y != 0 && fit_min_y < fit_rmse && fit_hit_y >= fit_range[0] && fit_hit_y <= fit_range[1] && is_have_target_det)
    {
        vld_fit_hit_y = fit_hit_y;
        vld_fit_slope_y = slope_y;
        vld_fit_intercept_y = Intercept_y;
        vector<float> temp_data = {(float)trigger_id, slope_y, Intercept_y, fit_min_y, fit_hit_y};
        // temp_data.insert(temp_data.end(), out_data.begin(), out_data.end());
        fec_fit_datay.push_back(temp_data);
    }
    else if (fit_min_y < fit_rmse && !is_have_target_det)
    {
        vector<float> temp_data = {(float)trigger_id, slope_y, Intercept_y};
        temp_data.insert(temp_data.end(), out_data_y.begin(), out_data_y.end());

        fec_fit_datay.push_back(temp_data);
    }
}

void fit_pro::det_data_fit(int num, float *fec_data_combine, float *fec_data_z, float *fec_data_combine_cluster, bool is_x_direction, float &vld_fit_hit, float &vld_fit_slope, float &vld_fit_intercept)
{
    vld_fit_hit = 0;
    vld_fit_slope = 0;
    vld_fit_intercept = 0;
    float fit_hit = 0;
    float slope = 0;
    float Intercept = 0;
    float rmse = 0;
    float fit_min = fit_rmse;
    double tan = 0;
    double fit_angle = 0;
    int data_size = det_layer_used;

    int k = 0;
    vector<float> out_data;
    for (int j = 0; j < num; j++)
    {
        float *current_fec_data_combine = fec_data_combine + k;
        float *current_fec_data_combine_cluster = fec_data_combine_cluster + k;
        float *current_fec_data_z;
        if (is_rotate)
        {
            current_fec_data_z = fec_data_z + k;
        }
        else
        {
            current_fec_data_z = fec_data_z;
        }

        k += data_size;
        VectorXd Y(data_size);
        MatrixXd X(data_size, 2);
        for (int n = 0; n < data_size; n++)
        {
            X(n, 0) = 1.0;
            X(n, 1) = current_fec_data_z[n] + z_alignment[n];
            if (is_x_direction)
            {

                Y(n) = current_fec_data_combine[n] + x_alignment[n];
            }
            else
            {
                Y(n) = current_fec_data_combine[n] + y_alignment[n];
            }
        }
        VectorXd beta = X.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Y);

        // // rmse
        rmse = 0;
        for (int z = 0; z < data_size; z++)
        {
            double fit_value = X.row(z) * beta; // 适用于多元
            rmse += (Y(z) - fit_value) * (Y(z) - fit_value);
        }
        rmse = std::sqrt(rmse / data_size);

        if (is_fit_control_angle)
        {
            tan = beta(1);
            fit_angle = std::atan(std::abs(tan)) * (180 / M_PI);
        }
        else
        {
            fit_angle = 0;
        }

        if (rmse < fit_rmse && rmse < fit_min && fit_angle < angle_deg)
        {
            out_data.clear();
            fit_min = rmse;
            // rmax = r2;
            slope = beta(1);
            Intercept = beta(0);
            for (int z = 0; z < data_size; z++)
            {
                out_data.push_back(Y(z));
            }

            for (int z = 0; z < data_size; z++)
            {
                out_data.push_back(X(z, 1));
            }
            if (is_write_event_json)
            {
                for (int z = 0; z < data_size; z++)
                {
                    out_data.push_back(current_fec_data_combine_cluster[z]);
                }
            }
            if (is_have_target_det)
                fit_hit = (fit_location + tz) * slope + Intercept;
        }
    }

    if (fit_hit != 0 && fit_min < fit_rmse && fit_hit >= fit_range[0] && fit_hit <= fit_range[1] && is_have_target_det)
    {
        vld_fit_hit = fit_hit;
        vld_fit_slope = slope;
        vld_fit_intercept = Intercept;
        vector<float> temp_data = {(float)trigger_id, slope, Intercept, fit_min, fit_hit};
        // temp_data.insert(temp_data.end(), out_data.begin(), out_data.end());
        if (is_x_direction)
        {
            fec_fit_datax.push_back(temp_data);
        }
        else
        {
            fec_fit_datay.push_back(temp_data);
        }
    }
    else if (fit_min < fit_rmse && !is_have_target_det)
    {
        vector<float> temp_data = {(float)trigger_id, slope, Intercept};
        temp_data.insert(temp_data.end(), out_data.begin(), out_data.end());

        if (is_x_direction)
        {
            fec_fit_datax.push_back(temp_data);
        }
        else
        {
            fec_fit_datay.push_back(temp_data);
        }
    }
}

void fit_pro::calculateSpatial(float *target_det_rz, float *target_det_chn, float *target_det_amp, float *target_det_clu, int target_det_size, float fit_hit, float fit_slope, float fit_intercept, float align_data, float &amp, float &chn, float &clu, float &x_err, float &y_err, bool is_x_direction)
{

    float min = eff_sigma;
    float chnTemp = 0;
    float ampTemp = 0;
    float clusterTemp = 0;
    float deta = 0;

    float *current_target_det_rz = target_det_rz;
    float *current_target_det_chn = target_det_chn;
    float *current_target_det_amp = target_det_amp;
    float *current_target_det_clu = target_det_clu;

    for (int j = 0; j < target_det_size; ++j)
    {
        deta = current_target_det_chn[j] - (fit_slope * (current_target_det_rz[j] + tz) + fit_intercept) + align_data;
        float abs_deta = abs(deta);
        if (abs_deta < abs(min))
        {
            min = deta;
            ampTemp = current_target_det_amp[j];
            clusterTemp = current_target_det_clu[j];
            chnTemp = current_target_det_chn[j];
        }
    }

    if (min != eff_sigma && chnTemp != 0 && ampTemp != 0)
    {
        if (is_x_direction)
        {
            rmsx.push_back(min);
            ampx.push_back(ampTemp);
            chnx.push_back(chnTemp);
            clusterx.push_back(clusterTemp);
            targe_vld_idx.push_back(trigger_id);
            x_err = min;
        }
        else
        {
            rmsy.push_back(min);
            ampy.push_back(ampTemp);
            chny.push_back(chnTemp);
            clustery.push_back(clusterTemp);
            targe_vld_idy.push_back(trigger_id);
            y_err = min;
        }
        amp = ampTemp;
        chn = chnTemp;
        clu = clusterTemp;
    }
}

void fit_pro::write_root_file()
{
    string filename = save_path + "adas_track_data.root";
    adas_track_file = new TFile(filename.c_str(), "recreate");
    adas_track_tree = new TTree("adas_track_data", "fit data");

    adas_track_tree->Branch("event_id", &event_id, "event_id/I");
    adas_track_tree->Branch("x_slope", &x_slope, "x_slope/F");
    adas_track_tree->Branch("x_Intercept", &x_Intercept, "x_Intercept/F");
    adas_track_tree->Branch("y_slope", &y_slope, "y_slope/F");
    adas_track_tree->Branch("y_Intercept", &y_Intercept, "y_Intercept/F");

    for (int i = 0; i < det_layer_used; i++)
    {
        string det_name = "detector" + to_string(i);
        adas_track_tree->Branch(det_name.c_str(), &detector_vld_data[i].x, "x/F:y/F:z/F");
    }

    cout << cBLUE << "root file write success" << cRESET << endl;
}

void fit_pro::write_target_root_file()
{
    string filename = save_path + "layer_" + to_string(target_det_idx + 1) + "_target_data.root";
    adas_target_file = new TFile(filename.c_str(), "recreate");
    adas_target_tree = new TTree("adas_target_data", "fit data");

    adas_target_tree->Branch("event_id", &event_id, "event_id/I");
    adas_target_tree->Branch("layer_num", &layer_num, "layer_num/I");
    adas_target_tree->Branch("det_idx", &det_idx, "det_idx/I");
    adas_target_tree->Branch("is_x_meet", &is_x_meet, "is_x_meet/I");
    adas_target_tree->Branch("is_y_meet", &is_y_meet, "is_y_meet/I");
    adas_target_tree->Branch("x_slope", &x_slope, "x_slope/F");
    adas_target_tree->Branch("x_Intercept", &x_Intercept, "x_Intercept/F");
    adas_target_tree->Branch("x_rmse", &x_rmse, "x_rmse/F");
    adas_target_tree->Branch("y_slope", &y_slope, "y_slope/F");
    adas_target_tree->Branch("y_Intercept", &y_Intercept, "y_Intercept/F");
    adas_target_tree->Branch("y_rmse", &y_rmse, "y_rmse/F");
    adas_target_tree->Branch("x_pos", &x_pos, "x_pos/F");
    adas_target_tree->Branch("y_pos", &y_pos, "y_pos/F");
    adas_target_tree->Branch("z_pos", &z_pos, "z_pos/F");
    adas_target_tree->Branch("x_pos_err", &x_pos_err, "x_pos_err/F");
    adas_target_tree->Branch("y_pos_err", &y_pos_err, "y_pos_err/F");
    adas_target_tree->Branch("x_amp", &x_amp, "x_amp/F");
    adas_target_tree->Branch("y_amp", &y_amp, "y_amp/F");
    adas_target_tree->Branch("x_strip_number", &x_strip_number, "x_strip_number/I");
    adas_target_tree->Branch("y_strip_number", &y_strip_number, "y_strip_number/I");

    cout << cBLUE << "target root file write success" << cRESET << endl;
}

void fit_pro::complete_adding()
{
    if (!is_have_target_det)
    {
        if (adas_track_tree != nullptr)
        {
            adas_track_file->cd();
            adas_track_tree->Write();
            adas_track_file->Flush();
            adas_track_file->Close();
        }
        else
        {
            cout << "adas_track_tree is not initialized!" << std::endl;
        }
    }
    else
    {
        if (adas_target_tree != nullptr)
        {
            adas_target_file->cd();
            adas_target_tree->Write();
            adas_target_file->Flush();
            adas_target_file->Close();
        }
        else
        {
            cout << "adas_target_tree is not initialized!" << std::endl;
        }
    }
}

void fit_pro::set_layer_used()
{
    for (int i = 0; i < valid_det_num; i++)
    {
        user_layer[i] = (use_set_layer_used >> i) & 1;
    }
    cout << "Layer used: " << endl;
    for (int i = 0; i < valid_det_num; i++)
    {
        if (user_layer[i])
        {
            cout << "Layer " << i << " ";
            cout << cGREEN << "used" << endl;
            cout << cRESET;
        }
    }
    cout << "========" << endl;
}

bool fit_pro::read_dec_file(string filename, vector<float> &dec)
{
    ifstream dec_file(filename);
    if (!dec_file.is_open())
    {
        cout << cRED << filename << " open failed" << endl;
        return false;
    }
    else
    {
        float tmp;
        while (dec_file >> tmp)
        {
            dec.push_back(tmp);
        }
        dec_file.close();
        return true;
    }
}

void fit_pro::read_alignment()
{
    vector<float> alignmentx, alignmenty, alignmentz;

    is_align = !(all_of(std::begin(alignment), std::end(alignment), [](float i)
                        { return i == 0; }));

    for (int i = 0; i < alignment.size(); i++)
    {
        if (i < valid_det_num)
        {
            alignmentx.push_back(alignment[i]);
        }
        else if (i < valid_det_num * 2 && i >= valid_det_num)
        {
            alignmenty.push_back(alignment[i]);
        }
        else if (i < valid_det_num * 3 && i >= valid_det_num * 2)
        {
            alignmentz.push_back(alignment[i]);
        }
    }

    tx = alignmentx[target_det_idx];
    ty = alignmenty[target_det_idx];
    tz = alignmentz[target_det_idx];

    for (int i = 0; i < alignmentx.size(); i++)
    {
        if (!(!user_layer[i] || (i == target_det_idx && is_have_target_det)))
        {
            x_alignment.push_back(alignmentx[i]);
            y_alignment.push_back(alignmenty[i]);
            z_alignment.push_back(alignmentz[i]);
        }
    }
    // cout << tx << " " << ty << " " << tz << endl;

    // for (int i = 0; i < x_alignment.size(); i++)
    // {
    //     cout << "detector " << i << " x alignment: " << x_alignment[i] << endl;
    //     cout << "detector " << i << " y alignment: " << y_alignment[i] << endl;
    //     cout << "detector " << i << " z alignment: " << z_alignment[i] << endl;
    // }
    cout << cGREEN << "Alignment ? " << is_align << endl;
    cout << cGREEN << "Alignment data read success!" << cRESET << endl;
    cout << "========" << endl;
}

void fit_pro::read_rotate()
{
    vector<Matrix<float, 3, 3>> rotate_matrix_wait;
    for (int i = 0; i < rotate_data.size(); i = i + 9)
    {
        Matrix<float, 3, 3> rotate_matrix_temp;
        rotate_matrix_temp << rotate_data[i], rotate_data[i + 1], rotate_data[i + 2],
            rotate_data[i + 3], rotate_data[i + 4], rotate_data[i + 5],
            rotate_data[i + 6], rotate_data[i + 7], rotate_data[i + 8];
        rotate_matrix_wait.push_back(rotate_matrix_temp);
    }

    rotate_matrix_target = rotate_matrix_wait[target_det_idx];

    for (int i = 0; i < rotate_matrix_wait.size(); i++)
    {
        if (!(!user_layer[i] || (i == target_det_idx && is_have_target_det)))
        {
            rotate_matrix.push_back(rotate_matrix_wait[i]);
        }
    }

    // for (int i = 0; i < rotate_matrix.size(); i++)
    // {
    //     cout << "detector " << i << " rotate matrix: " << endl;
    //     cout << rotate_matrix[i] << endl;
    // }

    // cout << rotate_matrix_target << endl;

    vector<bool> is_identity;
    for (const auto &matrix : rotate_matrix)
    {
        is_identity.push_back(is_identity_matrix(matrix));
    }

    is_rotate = !(all_of(std::begin(is_identity), std::end(is_identity), [](bool i)
                         { return i; }));

    cout << cGREEN << "Rotate ? " << is_rotate << endl;
    cout << cGREEN << "Rotate data read success!" << cRESET << endl;
    cout << "========" << endl;
}

bool fit_pro::is_identity_matrix(const Eigen::Matrix<float, 3, 3> &matrix)
{
    // 检查矩阵是否为方阵
    if (matrix.rows() != matrix.cols())
    {
        return false;
    }

    // 检查对角线元素是否全为1，非对角线元素是否全为0
    for (int i = 0; i < matrix.rows(); ++i)
    {
        for (int j = 0; j < matrix.cols(); ++j)
        {
            if (i == j)
            {
                if (matrix(i, j) != 1.0)
                {
                    return false;
                }
            }
            else
            {
                if (matrix(i, j) != 0.0)
                {
                    return false;
                }
            }
        }
    }

    return true;
}

void fit_pro::draw_region_amp()
{
    int nbins = 50; // 每个区域的bin数
    int max_index = 0;
    if (ADC_LSB == 1)
    {
        max_index = 10000;
    }
    else
    {
        max_index = 1000;
    }
    // for (int i = 0; i < region_number; ++i)
    // {
    //     std::string cname = "c_region_" + std::to_string(i + 1);
    //     std::string hname = "h_region_" + std::to_string(i + 1);
    //     std::string htitle = "Region " + std::to_string(i + 1) + " Amplitude;Amplitude;Counts";

    //     TCanvas *c = new TCanvas(cname.c_str(), htitle.c_str(), 800, 600);
    //     TH1D *h = new TH1D(hname.c_str(), htitle.c_str(), nbins, 0, max_index);

    //     // 填充直方图
    //     for (double amp : region_par_amp[i])
    //     {
    //         h->Fill(amp);
    //     }

    //     h->Draw();
    //     c->SaveAs((save_path + cname + ".png").c_str());
    // }

    int m = sqrt(region_number);

    for (int i = 0; i < m; i++)
    {
        std::string cname = "layer" + to_string(target_det_idx + 1) + "x_region_" + std::to_string(i + 1);
        std::string hname = "layer" + to_string(target_det_idx + 1) + "x_region_" + std::to_string(i + 1);
        std::string htitle = "x " + std::to_string(i + 1) + " Amplitude;Amplitude;Counts";

        TCanvas *c = new TCanvas(cname.c_str(), htitle.c_str(), 800, 600);
        TH1D *h = new TH1D(hname.c_str(), htitle.c_str(), nbins, 0, max_index);

        // 填充直方图
        for (double amp : target_region_effx[i])
        {
            h->Fill(amp);
        }

        h->Draw();
        c->SaveAs((save_path + cname + ".png").c_str());
    }

    for (int i = 0; i < m; i++)
    {
        std::string cname = "layer" + to_string(target_det_idx + 1) + "y_region_" + std::to_string(i + 1);
        std::string hname = "layer" + to_string(target_det_idx + 1) + "y_region_" + std::to_string(i + 1);
        std::string htitle = "y " + std::to_string(i + 1) + " Amplitude;Amplitude;Counts";

        TCanvas *c = new TCanvas(cname.c_str(), htitle.c_str(), 800, 600);
        TH1D *h = new TH1D(hname.c_str(), htitle.c_str(), nbins, 0, max_index);

        // 填充直方图
        for (double amp : target_region_effy[i])
        {
            h->Fill(amp);
        }

        h->Draw();
        c->SaveAs((save_path + cname + ".png").c_str());
    }
}

void fit_pro::draw_region_efficiency()
{
    vector<float> region_efficiency;
    vector<float> region_efficiency_x;
    vector<float> region_efficiency_y;
    for (int i = 0; i < sqrt(region_number); i++)
    {
        cout << cGREEN << "Region " << i + 1 << " efficiency: " << endl;
        outall << "Region" << i + 1 << "_efficiency ";
        if (is_region_eff)
        {
            float x = float(target_region_effx[i].size()) / float(ref_region_effx[i].size());
            float y = float(target_region_effy[i].size()) / float(ref_region_effy[i].size());
            region_efficiency_x.push_back(x);
            region_efficiency_y.push_back(y);
            cout << cGREEN << "X direction: " << x << endl;
            cout << cGREEN << "Y direction: " << y << endl;
            outall << "X: " << x << " Y: " << y << endl;
        }
        else
        {
            cout << cRED << "Region efficiency not calculated" << cRESET << endl;
        }
    }
    cout << endl;
    for (int i = 0; i < region_number; i++)
    {
        cout << cGREEN << "Region " << i + 1 << " efficiency: " << endl;
        outall << "Region" << i + 1 << "_efficiency ";
        if (is_region_eff)
        {
            float xy = float(target_region_par[i].size()) / float(ref_region_par[i].size());
            region_efficiency.push_back(xy);
            cout << cGREEN << xy << endl;
            outall << xy << endl;
        }
        else
        {
            cout << cRED << "Region efficiency not calculated" << cRESET << endl;
        }
    }

    int nrow = sqrt(region_number); // 行数
    int ncol = sqrt(region_number); // 列数
    double region = fit_range[1] - fit_range[0];
    double region_step = region / nrow; // 每个区域的步长

    // X方向效率热力图
    TCanvas *c1 = new TCanvas("c1", "X Direction Efficiency", 800, 700);
    TH2D *h2x = new TH2D("h2x", "X Direction Efficiency;X (mm);Efficiency", ncol, 0.5, (ncol + 0.5) * region_step, 1, 0, 1);
    for (int row = 0; row < 1; ++row)
    {
        for (int col = 0; col < ncol; ++col)
        {
            int idx = row * ncol + col;
            h2x->SetBinContent(col + 1, row + 1, region_efficiency_x[idx]);
        }
    }
    string xname = "layer" + to_string(target_det_idx + 1) + "_x_efficiency";
    h2x->SetMinimum(0);
    h2x->SetMaximum(1.0);
    gStyle->SetPaintTextFormat(".2f"); // 保留2位小数
    h2x->SetMarkerSize(2.0);           // 字体变大
    h2x->Draw("COLZ TEXT");
    h2x->SetTitle(xname.c_str());
    gStyle->SetOptStat(0); // 关闭统计框
    c1->SaveAs((save_path + xname + ".png").c_str());

    // Y方向效率热力图
    TCanvas *c2 = new TCanvas("c2", "Y Direction Efficiency", 800, 700);
    TH2D *h2y = new TH2D("h2y", "Y Direction Efficiency;Y (mm);Efficiency", ncol, 0.5, (ncol + 0.5) * region_step, 1, 0, 1);
    for (int row = 0; row < 1; ++row)
    {
        for (int col = 0; col < ncol; ++col)
        {
            int idx = row * ncol + col;
            h2y->SetBinContent(col + 1, row + 1, region_efficiency_y[idx]);
        }
    }
    string yname = "layer" + to_string(target_det_idx + 1) + "_y_efficiency";
    h2y->SetTitle(yname.c_str());
    h2y->SetMinimum(0);
    h2y->SetMaximum(1.0);
    gStyle->SetPaintTextFormat(".2f"); // 保留2位小数
    h2y->SetMarkerSize(2.0);           // 字体变大
    h2y->Draw("COLZ TEXT");
    gStyle->SetOptStat(0); // 关闭统计框
    c2->SaveAs((save_path + yname + ".png").c_str());

    // xy区域效率图
    TCanvas *c3 = new TCanvas("c1", "Region Efficiency Heatmap", 800, 700);
    TH2D *h2 = new TH2D("h2", "Region Efficiency;X (mm);Y (mm)", ncol, 0.5, (ncol + 0.5) * region_step, nrow, 0.5, (nrow + 0.5) * region_step);

    // 填充每个区域的效率
    for (int row = 0; row < nrow; ++row)
    {
        for (int col = 0; col < ncol; ++col)
        {
            int idx = row * ncol + col;
            h2->SetBinContent(col + 1, row + 1, region_efficiency[idx]); // y轴反向显示
        }
    }
    string name = "layer" + to_string(target_det_idx + 1) + "_region_efficiency";
    h2->SetMinimum(0);
    h2->SetMaximum(1.0);
    gStyle->SetPaintTextFormat(".2f"); // 保留2位小数
    h2->SetMarkerSize(2.0);            // 字体变大
    h2->Draw("COLZ TEXT");             // COLZ: 颜色, TEXT: 显示数值
    h2->SetTitle(name.c_str());
    gStyle->SetOptStat(0); // 关闭统计框
    c3->SaveAs((save_path + name + ".png").c_str());
}

void fit_pro::draw_hit_strip(vector<float> ss, bool xory)
{
    string name = "hit_strip";
    if (xory)
    {
        name = "hit_strip_layer_x" + to_string(target_det_idx + 1);
    }
    else
    {
        name = "hit_strip_layer_y" + to_string(target_det_idx + 1);
    }

    TCanvas *c1 = new TCanvas("c1", "Hit Strip", 800, 600);
    TH1F *h1 = new TH1F("h1", "Hit Strip", 120, 0, 600);
    h1->SetXTitle(xory ? "X (mm)" : "Y (mm)");
    h1->SetYTitle("Counts");
    h1->SetTitle(name.c_str());
    gStyle->SetOptStat(0); // 关闭统计框

    for (const auto &ss_vec : ss)
    {
        h1->Fill(ss_vec);
    }

    h1->Draw();

    c1->SaveAs((save_path + name + ".png").c_str());
}

void fit_pro::draw_hit_range_xy()
{
    int N = x_all.size();

    double x[N], y[N];

    for (int i = 0; i < N; ++i)
    {
        x[i] = x_all[i];
        y[i] = y_all[i];
    }

    TCanvas *c1 = new TCanvas("c1", "Hit Map", 800, 600);
    TGraph *gr = new TGraph(N, x, y);
    gr->SetTitle("2D Hit Map;X (mm);Y (mm)");
    gr->SetMarkerStyle(20);
    gr->SetMarkerSize(0.7);
    gr->Draw("AP"); // A:坐标轴, P:点

    string name = "hit_map_layer" + to_string(target_det_idx + 1);
    c1->SaveAs((save_path + name + ".png").c_str());
}

void fit_pro::draw_spatial(vector<float> rms, bool xory)
{
    double max_index, min_index;
    double fit_index;
    if (is_align)
    {
        max_index = 5;
        min_index = -5;
        fit_index = 3;
    }
    else
    {
        max_index = 20;
        min_index = -20;
    }

    TCanvas *c2 = new TCanvas("c2", "spitial_resolution", 800, 600);
    string name = (xory ? "x" : "y");
    name = name + "_spatial_layer" + to_string(target_det_idx + 1);
    TH1F *spatial = new TH1F(name.c_str(), "Histogram from rms", 100, min_index, max_index);

    for (const auto &rms_vec : rms)
    {
        spatial->Fill(rms_vec);
    }

    TF1 *fitFcn = new TF1("fitFcn", "gaus(0) + gaus(3)", -fit_index, fit_index);
    fitFcn->SetNpx(1600);

    // Set initial parameter values
    fitFcn->SetParameters(spatial->GetMaximum(), spatial->GetMean(), spatial->GetRMS(),
                          spatial->GetMaximum() / 2, spatial->GetMean(), spatial->GetRMS() * 2);

    // Fit histogram
    string xais = (xory ? "#Deltax (mm)" : "#Deltay (mm)");
    spatial->Fit("fitFcn", "R");
    spatial->GetXaxis()->SetTitle(xais.c_str());
    spatial->GetYaxis()->SetTitle("Counts");
    spatial->GetXaxis()->CenterTitle();
    spatial->GetYaxis()->CenterTitle();

    double lowerBound = -fit_index;
    double upperBound = fit_index;

    while (fitFcn->GetParameter(0) < 0 || fitFcn->GetParameter(3) < 0)
    {
        lowerBound += 0.5;
        upperBound -= 0.5;
        delete fitFcn;
        fitFcn = new TF1("fitFcn", "gaus(0) + gaus(3)", lowerBound, upperBound);
        fitFcn->SetNpx(1600);
        fitFcn->SetParameters(spatial->GetMaximum(), spatial->GetMean(), spatial->GetRMS(),
                              spatial->GetMaximum() / 2, spatial->GetMean(), spatial->GetRMS() * 2);
        spatial->Fit("fitFcn", "R");
    }

    spatial->Draw();

    gStyle->SetOptFit(1111); // 打开拟合参数信息
    c2->Update();            // 更新画布以确保统计框已创建

    TF1 *gaus1 = new TF1("gaus1", "gaus", -fit_index, fit_index);
    TF1 *gaus2 = new TF1("gaus2", "gaus", -fit_index, fit_index);

    // Set their parameters
    gaus1->SetParameters(fitFcn->GetParameter(0), fitFcn->GetParameter(1), fitFcn->GetParameter(2));
    gaus2->SetParameters(fitFcn->GetParameter(3), fitFcn->GetParameter(4), fitFcn->GetParameter(5));
    gaus1->SetLineColor(kBlue);
    int brown = TColor::GetColor("#8B4513");
    gaus2->SetLineColor(brown);
    // Draw them on the same canvas
    gaus1->Draw("same");
    gaus2->Draw("same");

    // 获取两个高斯函数的参数
    double A1 = fitFcn->GetParameter(0);     // 高斯1的幅度
    double mean1 = fitFcn->GetParameter(1);  // 高斯1的均值
    double sigma1 = fitFcn->GetParameter(2); // 高斯1的标准差

    double A2 = fitFcn->GetParameter(3);     // 高斯2的幅度
    double mean2 = fitFcn->GetParameter(4);  // 高斯2的均值
    double sigma2 = fitFcn->GetParameter(5); // 高斯2的标准差

    // 计算core sigma
    double eff_sigma = A1 > A2 ? sigma1 : sigma2;

    A1 = A1 * abs(sigma1);
    A2 = A2 * abs(sigma2);

    // 计算加权均值
    double weighted_mean = (A1 * abs(mean1) + A2 * abs(mean2)) / (A1 + A2);

    // 计算加权 RMS
    double weighted_rms = std::sqrt(
        (A1 * (sigma1 * sigma1 + mean1 * mean1) + A2 * (sigma2 * sigma2 + mean2 * mean2)) / (A1 + A2) - weighted_mean * weighted_mean);

    cout << cMAGENTA << "Sigma of Gaussian: " << eff_sigma << endl;
    cout << cMAGENTA << "Weighted RMS: " << weighted_rms << endl;

    // Save canvas as PDF file
    string filename;
    string sigmaname;
    int angle_set = angle_deg;

    if (xory)
    {
        filename = save_path + to_string(target_det_idx + 1) + "_spatial_x" + (is_fit_control_angle ? "_angle" + to_string(angle_set) : "") + (is_align ? "_align" : "") + ".png";
        sigmaname = "Xspatial_" + to_string(target_det_idx + 1);
    }
    else
    {
        filename = save_path + to_string(target_det_idx + 1) + "_spatial_y" + (is_fit_control_angle ? "_angle" + to_string(angle_set) : "") + (is_align ? "_align" : "") + ".png";
        sigmaname = "Yspatial_" + to_string(target_det_idx + 1);
    }

    outall << sigmaname << " " << eff_sigma << endl;
    outall << sigmaname << " weighted_rms: " << weighted_rms << endl;
    cout << cRED << "detector " << target_det_idx + 1 << " " << (xory ? "x" : "y") << " direction eff_sigma: " << eff_sigma << cRESET << endl;
    c2->SaveAs(filename.c_str());

    spatial->Reset();
    delete spatial;
    delete fitFcn;
    delete gaus1;
    delete gaus2;
    delete c2;
}

void fit_pro::draw_energy(vector<float> energy, bool xory)
{
    int max_index = 0;
    if (ADC_LSB == 1)
    {
        max_index = 10000;
    }
    else
    {
        max_index = 1000;
    }
    string name = (xory ? "x" : "y");
    name = name + "_energy_layer" + to_string(target_det_idx + 1);
    TH1F *muonenergy = new TH1F(name.c_str(), "Histogram from muon deposit energy", 100, -5, max_index);

    for (const auto &energy_vec : energy)
    {
        muonenergy->Fill(energy_vec);
    }

    TF1 *fitFcn = new TF1("fitFcn", "[0]*TMath::Landau(x,[1],[2])*TMath::Gaus(x,[3],[4])", -5, max_index);
    fitFcn->SetNpx(500);

    // Fit histogram
    string xais = ((ADC_LSB == 1) ? "ADC code" : "Charge (fC)");
    muonenergy->Fit("fitFcn", "R");
    muonenergy->GetXaxis()->SetTitle(xais.c_str());
    muonenergy->GetYaxis()->SetTitle("Counts");
    muonenergy->GetXaxis()->CenterTitle();
    muonenergy->GetYaxis()->CenterTitle();

    // 设置初始参数值
    fitFcn->SetParameter(0, muonenergy->GetMaximum()); // 高度
    fitFcn->SetParameter(1, muonenergy->GetMean());    // 朗道的中心
    fitFcn->SetParameter(2, muonenergy->GetRMS());     // 朗道的宽度
    fitFcn->SetParameter(3, muonenergy->GetMean());    // 高斯的中心
    fitFcn->SetParameter(4, muonenergy->GetRMS());     // 高斯的宽度

    // 拟合直方图
    muonenergy->Fit("fitFcn", "R");

    // 创建一个新的画布
    TCanvas *c2 = new TCanvas("c2", "Canvas", 800, 600);

    // 在新的画布上绘制原始直方图和拟合函数
    muonenergy->Draw();
    fitFcn->Draw("same"); // 在原始直方图上添加拟合函数的图像

    double energy_mip = fitFcn->GetParameter(1);
    // 保存画布为PDF文件
    // Save canvas as PDF file
    string filename;
    string energyname;
    int angle_set = angle_deg;

    if (xory)
    {
        filename = save_path + to_string(target_det_idx + 1) + "_energy_x" + (is_fit_control_angle ? "_angle" + to_string(angle_set) : "") + (is_align ? "_align" : "") + ".png";
        energyname = "Xmuonenergy_" + to_string(target_det_idx + 1);
    }
    else
    {
        filename = save_path + to_string(target_det_idx + 1) + "_energy_y" + (is_fit_control_angle ? "_angle" + to_string(angle_set) : "") + (is_align ? "_align" : "") + ".png";
        energyname = "Ymuonenergy_" + to_string(target_det_idx + 1);
    }

    outall << energyname << " " << energy_mip << endl;
    cout << cRED << "detector " << target_det_idx + 1 << " " << (xory ? "x" : "y") << " direction mip: " << energy_mip << cRESET << endl;
    c2->SaveAs(filename.c_str());

    muonenergy->Reset();
    delete muonenergy;
    delete fitFcn;
    delete c2;
}

void fit_pro::draw_all_energy(vector<float> energy)
{

    int max_index = 0;
    if (ADC_LSB == 1)
    {
        max_index = 10000;
    }
    else
    {
        max_index = 1000;
    }
    string name = "xy_energy_layer" + to_string(target_det_idx + 1);
    TH1F *muonenergy = new TH1F(name.c_str(), "Histogram from muon deposit energy", 100, -5, max_index);

    for (const auto &energy_vec : energy)
    {
        muonenergy->Fill(energy_vec);
    }

    TF1 *fitFcn = new TF1("fitFcn", "[0]*TMath::Landau(x,[1],[2])*TMath::Gaus(x,[3],[4])", -5, max_index);
    fitFcn->SetNpx(500);

    // 设置初始参数值
    fitFcn->SetParameter(0, muonenergy->GetMaximum()); // 高度
    fitFcn->SetParameter(1, muonenergy->GetMean());    // 朗道的中心
    fitFcn->SetParameter(2, muonenergy->GetRMS());     // 朗道的宽度
    fitFcn->SetParameter(3, muonenergy->GetMean());    // 高斯的中心
    fitFcn->SetParameter(4, muonenergy->GetRMS());     // 高斯的宽度

    // 拟合直方图
    muonenergy->Fit("fitFcn", "R");

    // Fit histogram
    string xais = ((ADC_LSB == 1) ? "ADC code" : "Charge (fC)");
    muonenergy->Fit("fitFcn", "R");
    muonenergy->GetXaxis()->SetTitle(xais.c_str());
    muonenergy->GetYaxis()->SetTitle("Counts");
    muonenergy->GetXaxis()->CenterTitle();
    muonenergy->GetYaxis()->CenterTitle();

    // 创建一个新的画布
    TCanvas *c2 = new TCanvas("c2", "Canvas", 800, 600);

    // 在新的画布上绘制原始直方图和拟合函数
    muonenergy->Draw();
    fitFcn->Draw("same"); // 在原始直方图上添加拟合函数的图像

    double energy_mip = fitFcn->GetParameter(1);
    cout << cMAGENTA << "XYmip of land-Gaussian: " << energy_mip << endl;

    // 保存画布为PDF文件
    string filename = save_path + to_string(target_det_idx + 1) + (is_align ? "_align" : "") + "_energy_xy.png";
    string energyname = "XYmuonenergy_" + to_string(target_det_idx + 1);
    outall << energyname << " " << energy_mip << endl;
    c2->SaveAs(filename.c_str());

    muonenergy->Reset();
    delete muonenergy;
    delete fitFcn;
    delete c2;
}

void fit_pro::draw_cluster(vector<float> cluster, bool xory)
{
    string name = (xory ? "x" : "y");
    name = name + "_hit_strip_num_layer" + to_string(target_det_idx + 1);
    TH1F *stripnum = new TH1F(name.c_str(), "Histogram from muon hit strip  number", 30, 0, 30);

    for (const auto &cluster_vec : cluster)
    {
        stripnum->Fill(cluster_vec);
    }

    // Create a new canvas
    TCanvas *c1 = new TCanvas("c1", "Canvas", 800, 600);

    // Draw the histogram on the canvas
    stripnum->Draw();

    // Save the canvas as a PNG file
    string filename;
    if (xory)
    {
        filename = save_path + to_string(target_det_idx + 1) + "_stripnum_x.png";
    }
    else
    {
        filename = save_path + to_string(target_det_idx + 1) + "_stripnum_y.png";
    }
    c1->SaveAs(filename.c_str());

    // Clean up
    stripnum->Reset();
    delete stripnum;
    delete c1;
}

void fit_pro::out_result()
{
    string name;
    int angle_set = angle_deg;
    name = "detector_" + to_string(target_det_idx + 1) + "_";

    if (is_align)
    {
        name = name + "align_";
    }
    else
    {
        name = name;
    }

    if (is_fit_control_angle)
    {
        name = name + "angle_deg" + to_string(angle_set) + "_";
    }
    else
    {
        name = name;
    }

    outall.open(save_path + name + "result.txt", ios::app);
    if (!outall.is_open())
    {
        cout << cRED << "Failed to open file: " << (save_path + name + "result.txt") << endl;
    }
}

void fit_pro::write_event_json_file(int outevt)
{
    event_json.clear();
    event_json["type"] = "event";
    event_json["id"] = 35;
    event_json["event_id"] = trig_event;
    event_json["tracks"][0]["start_point"] = {
        kx * 600 + bx, ky * 600 + by, 600};
    event_json["tracks"][1]["start_point"] = {
        kx * -100 + bx, ky * -100 + by, -100};

    TVector3 vec(kx, ky, 1);
    vec.SetMag(1);
    event_json["tracks"][0]["direction"] = {
        vec.x(), vec.y(), vec.z()};
    event_json["tracks"][0]["angle_deg"] = {
        round(vec.Theta() / 3.1415 * 180 * 100) / 100, round(vec.Phi() / 3.1415 * 180 * 100) / 100};

    TVector3 vec1(kx, ky, 1);
    vec1.SetMag(1);
    event_json["tracks"][1]["direction"] = {
        vec1.x(), vec1.y(), vec1.z()};
    event_json["tracks"][1]["angle_deg"] = {
        round(vec1.Theta() / 3.1415 * 180 * 100) / 100, round(vec1.Phi() / 3.1415 * 180 * 100) / 100};

    int k = 0;
    for (int i = det_layer_used - 1; i >= 0; i--)
    {
        event_json["hits"][k] = {
            kx * z[i] + bx, ky * z[i] + by};
        int numx = x_strip_num[i] / 2;
        int numy = y_strip_num[i] / 2;
        for (int j = 0; j < x_strip_num[i]; j++)
        {
          event_json["x_cells"][k][j] = int(x_strip[i] / cSTRIP_WIDTH - numx) + j;
        }
        for (int j = 0; j < y_strip_num[i]; j++)
        {
          event_json["y_cells"][k][j] = int(y_strip[i] / cSTRIP_WIDTH - numy) + j;
        }
        k++;
    }

    string js_filename = save_path_event + "/event" + "/event_" + to_string(outevt) + ".json";
    ofstream js_file(js_filename);
    js_file << event_json << endl;
    js_file.close();
}

void fit_pro::write_status_json_file()
{
    status_json.clear();
    status_json["type"] = "event";
    status_json["id"] = 35;
    status_json["daq"] = "on";
    status_json["event_rate"] = 21.6;
    for (int i = 0; i < 4; i++)
    {
        status_json["hv_supply"][2 * i]["name"] = "hv" + to_string(i + 1) + "_cathode";
        status_json["hv_supply"][2 * i]["voltage"] = 720;
        status_json["hv_supply"][2 * i]["current"] = 0.001;
        status_json["hv_supply"][2 * i + 1]["name"] = "hv" + to_string(i + 1) + "_mesh";
        status_json["hv_supply"][2 * i + 1]["voltage"] = 560;
        status_json["hv_supply"][2 * i + 1]["current"] = 0.001;
    }
    string js_filename = save_path_event + "/status/" + "status_0.json";
    ofstream js_file(js_filename);
    js_file << status_json << endl;
    js_file.close();
}
